* Define program to estimate the damages
* This program is used inside the next program "ricardian_program"
capture program drop damages_precip
program define damages_precip, rclass

	syntax [, deficit(string)  edd(string) ///
			knots_d(real 0) knots_g(real 0) knots_e(real 0) ///
			rent_hat(string)]
	qui {
	local knots_`deficit' `knots_d'
	*local knots_`surplus' `knots_s'
	local knots_gdd `knots_g'
	local knots_`edd' `knots_e'

	* Use same knot locations as in my regression
	foreach var of varlist `deficit' gdd `edd' {
	local knot1=`var'_knots[1,1]
	local knot2=`var'_knots[1,2]
	local knot3=`var'_knots[1,3]
	local knot4=`var'_knots[1,4]
	local knot5=`var'_knots[1,5]
	local knot6=`var'_knots[1,6]
	capture drop `var'_spline*
		if `knots_`var''==3 {
		mkspline `var'_spline = `var', cubic knots(`knot1'  `knot2'  `knot3')
		}
		if `knots_`var''==4 {
		mkspline `var'_spline = `var', cubic knots(`knot1'  `knot2'  `knot3'  `knot4')
		}
		if `knots_`var''==5 {
		mkspline `var'_spline = `var', cubic knots(`knot1'  `knot2'  `knot3'  `knot4' `knot5')
		}
		if `knots_`var''==6 {
		mkspline `var'_spline = `var', cubic knots(`knot1'  `knot2'  `knot3'  `knot4' `knot5' `knot6')
		}
	}
	
	estimates restore precip_results
	predict ln_rent_hat_new, xb
	gen rent_hat_new=exp(ln_rent_hat_new+e(rmse)^2/2)
	
	gen chg_rent=rent_hat_new - `rent_hat'
	gen ln_chg_rent= ln_rent_hat_new - ln_rent_hat
	gen perc_chg_rent=chg_rent/`rent_hat'
	
	* Total Damages
	gen damages=chg_rent*croplandNon_acres
	total damages
	return scalar total_damages=_b[damages]
	} // end of qui

end


capture program drop ricardian_program_precip
program ricardian_program_precip, eclass
tempname parameters 

	syntax [, deficit(string) edd(string) ///
			knots_d(real 0) knots_g(real 0) knots_e(real 0) ///
			reg_year(real 0) soils_controls(string) scenario(string) model(string) wild(real 0)]

tempvar rent_hat total_rent

if `wild'==1{
capture drop temp pos wildresid wildy
bys stfips_clust: gen temp = uniform() 
by stfips_clust: gen pos = (temp[1] < .5) 
gen wildresid = epshat * (2*pos - 1) 
gen wildy = xb + wildresid 
}		

capture drop `deficit'_spline*
mkspline `deficit'_spline = `deficit', cubic nknots(`knots_d') displayknots
matrix `deficit'_knots=r(knots)

capture drop gdd_spline*
mkspline gdd_spline = gdd, cubic nknots(`knots_g') displayknots
matrix gdd_knots=r(knots)

capture drop `edd'_spline*
mkspline `edd'_spline = `edd', cubic nknots(`knots_e') displayknots
matrix `edd'_knots=r(knots)		

if `wild'==0{			
reg ln_cash_rent `deficit'_spline* gdd_spline* `edd'_spline* `soils_controls' i.stfips ///
			[aw=croplandNon_acres], vce(cluster stfips_clust)
est sto precip_results
* Predicted value and residual used for wild bootstrap
predict xb if e(sample), xb
predict epshat if e(sample), resid
}

if `wild'==1{			
reg wildy `deficit'_spline* gdd_spline* `edd'_spline* `soils_controls' i.stfips ///
			[aw=croplandNon_acres], vce(cluster stfips_clust)
est sto precip_results
}

qui {
*-------------------------------------------------------------------------------
* Estimate Climate Change Damages
*-------------------------------------------------------------------------------
predict ln_rent_hat, xb
gen `rent_hat'=exp(ln_rent_hat+e(rmse)^2/2)
drop ln_rent_hat

* total rent in the region
gen `total_rent'=`rent_hat'*croplandNon_acres
total `total_rent'
scalar base_rent=_b[`total_rent']
scalar ln_base_rent=ln(base_rent)


*******************************
* Total Damages
*******************************
preserve
	foreach var of varlist `deficit' gdd `edd' {
		summ `var' chg_`var'_rcp45AVG
		replace `var'=`var' + chg_`var'_`scenario'`model'
	}
	damages_precip, deficit("`deficit'") edd("`edd'") ///
			knots_d(`knots_d') knots_g(`knots_g') knots_e(`knots_e') ///
			rent_hat(`rent_hat')
	return list, all
	
	if `wild'==0{	
	save "..\temp\precip_damages_total_`scenario'`model'", replace
	}

restore

* Relative change in rent
local perc_chg_rent_total=r(total_damages)/base_rent
display "`perc_chg_rent_total'"
* Change in log rent
local ln_new_rent=r(total_damages)+base_rent
local chg_ln_rent_total= ln(`ln_new_rent') - ln_base_rent
display "`chg_ln_rent_total'"
* Relative change in rent implied by log change
display exp(`chg_ln_rent_total') - 1

*------------------------------------------------------------------------------------
/* If I calculate the percent change in rent due to climate damages
   for each variable separately, then the perecent changes do not all
   add to 1 perfectly due to the log tranformation. But note that the change
   in natural log only approximates a percent change and this approximation
   is less accurate for larger changes. Above I have shown two alternative 
   ways of calculating the percent change in rent and get the same answer either way.
   To decompose the total change in rent and find the share of damages due to each
   source, I calculate the share of the total change in log rent because these 
   shares add to 1. */
*------------------------------------------------------------------------------------

*******************************
* Only Precip
*******************************
preserve
	foreach var of varlist `deficit'  {
		replace `var'=`var' + chg_`var'_`scenario'`model'
	}
	damages_precip, deficit("`deficit'") edd("`edd'") ///
			knots_d(`knots_d') knots_g(`knots_g') knots_e(`knots_e') ///
			rent_hat(`rent_hat')
	return list, all
	
	if `wild'==0{	
	save "..\temp\precip_damages_precip_`scenario'`model'", replace
	}
	
restore

* Relative change in rent
local perc_chg_rent_deficit=r(total_damages)/base_rent
* Proportion due to deficit
local share_deficit = (ln(r(total_damages)+base_rent) - ln_base_rent) /`chg_ln_rent_total'


*******************************
* Only Heat Stress
*******************************
preserve
	foreach var of varlist  gdd `edd' {
		replace `var'=`var' + chg_`var'_`scenario'`model'
	}
	damages_precip, deficit("`deficit'") edd("`edd'") ///
			knots_d(`knots_d') knots_g(`knots_g') knots_e(`knots_e') ///
			rent_hat(`rent_hat')
	return list, all

	if `wild'==0{	
	save "..\temp\precip_damages_heat_`scenario'`model'", replace
	}

restore

* Relative change in rent
local perc_chg_rent_heat=r(total_damages)/base_rent
* Proportion due to heat
local share_heat = (ln(r(total_damages)+base_rent) - ln_base_rent) /`chg_ln_rent_total'


} // end of qui

matrix `parameters'= (`perc_chg_rent_total',`perc_chg_rent_deficit',`perc_chg_rent_heat', ///
					`share_deficit',`share_heat')

matrix colnames `parameters' = rel_total_chg rel_precip_chg rel_heat_chg ///
				share_precip share_heat 					
					
ereturn post `parameters'
end
